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ABSTRACT 

We present a new measurement of the mass of the MiUcy Way (MW) based on observed properties of its 
largest satellite galaxies, the Magellanic Clouds (MCs), and an assumed prior of a ACDM universe. The 
large, high-resolution Bolshoi cosmological simulation of this universe provides a means to statistically sample 
the dynamical properties of bright satellite galaxies in a large population of dark matter halos. The observed 
properties of the MCs, including their circular velocity, distance from the center of the MW, and velocity within 
the MW halo, are used to evaluate the likelihood that a given halo would have each or all of these properties; the 
posterior PDF for any property of the MW system can thus be constructed. This method provides a constraint 
on the MW virial mass, L2^Q4(stat.)^jj 3(sys.) x 10'~Mq (68% confidence), which is consistent with recent 
determinations that involve very different assumptions. In addition, we calculate the posterior PDF for the 
density profile of the MW and its satellite accretion history. Although typical satellites of 1O'~M0 halos are 
accreted over a wide range of epochs over the last 10 Gyr, we find a ^72% probability that the Magellanic 
Clouds were accreted within the last Gyr, and a 50%probability that they were accreted together 
Subject headings: Galaxy: formation, fundamental parameters, halo — galaxies: dwarf, Magellanic Clouds, 
evolution — dark matter 



1. INTRODUCTION 

The contents of the Milky Way (MW) Galaxy and the 
satellites that fall under its dynamical spell provide a unique 
testbed for theories of galaxy formation and cosmology. De- 
tailed observations of resolved stars, including proper mo- 
tions, allow the mass distribution of the galaxy to be measured 
with ever higher precision (e.g. iKallivavalil et alj I2006ct : 
iPiateket al]l2008h . Satellite galaxies within the MW have 
been detected with luminosities thr ee orders of magnitud e 
smaller than in external galaxies (e.g. lBelokurov et alj|2007l) . 
In addition, determining the detailed phase space distribu- 
tion of the MW galaxy is critical for interpreting the results 
of experim ents to directly or indirectly detect particle dark 
matter re.g.lStrig ari & Trotta"2009':'Vog elsbergereTan i2009l: 
iLisanti et al.ll201CkiKuhlen et al.. 2010) . A full understanding 
of the MW's place in the Universe requires not only detailed 
knowledge of its mass distribution and formation history, but 
also a sense of how this one well-studied system fits into the 
full cosmological context. 

A variety of methods have been used to put limits 
on the MW mass, ra nging from stellar dynamics and 

J2002I; 'B attagUa et al] 
j2p07; Xue et alJl2008F 
20IOI) to the dynamics 
20081) . For example. 



dynamics of satellites ( Klvpin et al 



20051: iDehnen et al.l l 2006l: ISmith et al 
Watkins et alj 120 1 oi 'Gnedin et al 



of the local group (Li & White 
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iBattagUa etall (12005b and IXue et alj (120081) performed a 
Jeans analysis of measurements of the radial velocity disper- 
sion profile from satellite galaxies, globular clusters, and blue 
horizontal b ranch halo s t ars to estimate the MW radial den- 
sity profile. ISmith et akl ( 120071) used measurements of high- 
velocity halo stars to estimate the MW escape speed assuming 
an NEW profile. iLi & Whit3 ( 120081) took a complementary 
approach, using the relative position and velocity of the MW 
and M3 1 , along with the age of the universe, to infer prop- 
erties of the obits of the MW-M31 system, which provides 
constraints on the total mass. 

Similarly, there have been a range of studies on the dy- 
namical state of the Magellanic Clouds (MCs). Until re- 
cently, the standard picture was that the M Cs were o bjects that 
have been orbiting the MW for some time dMurai & Fujimot^ 
[19801 iGardiner et al.|[T99l . This picture was motivated in 
part by the presence of the Magellanic Stream, a filament 
of gas extending 150° across the sky. Because it is appar- 
ently spatially and chemically associated with the Magellanic 
Clouds, it has often been interpreted as a tidal tail and taken 
as an indication that the satellites have been around f or sev- 
eral G yr (see, e.g., Lin & Lvnden-Bell 1982; Connor s et alj 
120061) . This picture has recently come under fire, largely as a 
result of detailed measurements of the three-dimensional ve- 
locity of the MCs: they are observed to have high velocities 
not aligned with the Magellanic Stream, indicating that they 
are not in virial equilibrium (and suggesting alte rnative for- 
mation methods for the Magellanic Stream; see Bes la et al] 
[2007. ,2010,) . Similarly, there is a growing consensus that the 
MCs accreted as a group: evid ence for this comes from both 
their proximity in phase space ( IKallivavalil et alj|2006al) . and 
the resu h that simulated subhalo s in general tend to accrete in 
groups (lD'Onghia&Lakdl2008h . 

In this work, we take a new approach to measure the 
mass and assembly of the MW. N-body simulations of 
dark matter structures in a ACDM universe have been 
very successful at reproducing the observed clustering 
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of galaxies (e.g. iKravtsov et alj I2004t iConrov et all 120061: 
rrruiillo-G omez et al.ll2010lK In two companion papers, we 
show tiiat the full probability distribution (PDF) for the num- 
ber of bright satellites around MW-lumino sity hosts predicted 
by high resolution numerical simulations jBusha et al.ll2010l) 
is in excel lent agree ment with measur ements from the SDSS 
jLiu et all [2010; To llerud et al.ll2011h . This provides evi- 
dence that such cosmological simulations realistically repre- 
sent galaxy halos and their satellites, and that they sample the 
underlying probability distribution for the properties of ha- 
los and subhalos in our Universe. These simulated halo cat- 
alogs therefore constitute a highly informative prior PDF for 
the parameters of any particular galaxy system. In this work 
we show that combining this prior with basic data about the 
two most massive MW satellites — their masses, velocities, 
and positions — provides interesting constraints on the MW 
mass, the distribution of mass within the MW system, and the 
system's assembly history. Although we find that the sim- 
ulation used here provides only a relatively sparse sampling 
of the underlying PDF, this is the first time that the statistics 
to study the MW system in this way have been available at 
all. As high-resolution cosmological simulations probe ever 
larger volumes, this approach will have increasing power and 
applicability. 

2. SIMULATIONS 

Statistical inference from halo dynamical histories requires 
a large, unbiased set of dark matter halos which samples the 
full range of cosmologically appropriate formati on scenarios. 
Here, we use halos from the Bolshoi simulation jKlvpin et al.l 
I20T0I) . which modeled a 250 /r'Mpc comoving box with 
n,n = 0.27, = 0.73, erg = 0.82, « = 0.95, and h = 0.7. 
The simulation volume contained 2048^ particles, each with 
mass 1.15 X lQ^h~^M (T,, and was run using the ART code 
( IKravtsov et al.l [19971) . Hal os and subhalos were iden tified 
using the BDM al gorithm jKlvpin & Holtzmaril II997I) : see 
iKlvpin et al.l ( 120 1 Oh for details. One unique aspect of this sim- 
ulation is the high spatial resolution, which is resolved down 
to a physical scale of 1 /2"'kpc. This improves the tracking 
of halos as they merge with and are disrupted by larger ob- 
jects, allowing them to be followed even as they pass near the 
core of the halo. The resulting halo catalog is nearly com- 
plete for objects down to a circular velocity of v^ax = 50 km 
i"'. While the overall simulation suffers from incomplete- 
ness in satellites at the ~ 20% level at 50km i"' , the incom- 
pleteness is strongly dependent on host mass. Host halos with 
Mvir = 0.5-3 X 1O'~/2~'M0 seem to be missing fewer than 
^ 10% of their subhalos. Because we are primarily interested 
in these lower mass objects, the amount of incompleteness is 
small and can be ignored. 

The large volume probed results in a sample of 2. 1 million 
simulated galaxy halos at the present epoch, including more 
than 100,000 halos massive enough to host at least one re- 
solved subhalo. We can increase this number further by con- 
sidering halos identified at different epochs to be independent 
objects representative of local systems: we use halos from 60 
simulation snapshots out to redshift 0.25. Throughout, we de- 
fine "hosts" as halos that are not within the virial radius of 
a larger halo, and "satellites" as any object within 300 kpc 
of a host. This value is chosen to be roughly half the dis- 
tance between the MW and M3 1 . Note that our final results 
are largely independent of this choice, because we will later 
impose more stringent requirements on the distance from the 



satellites to the center of their host halo. 

3. OBSERVATIONS 

We now consider the massive subhalo population of the 
MW that is modeled by the Bolshoi simulation: objects with 
Vmax > 50 km i"'. The two brightest MW satellite galax- 
ies, the LMC and SMC, have both been measured to have 
maximum circular velocities v^ax ^ 60 km s~^ w ith magni- 
tudes My = -18.5 and -17 . 1, respectivelv ([van der M arel et al] 
I2OOI IStanimirovicetal] I2004I: Ivan den Berghn200(]|) . The 
next brightest satellite is Sagittarius, some 4 magnitudes dim- 
mer, with Vmax ^ 20 km ^"^( iStrigari et alj|2007l) . Similar con- 
straints, albeit with larger error bars, can be made for the 
other bright classical satellites. The census of nearby objects 
brighter than My ^ -8 should be complete well beyond the 
MW virial radius dWalsh et alj|2009l: iTollerud et alj|2008l) . It 
is therefore a robust statement that the MW has exactly two 
satellites with v^ax ^ 50km . 

Applying this selection criterion to the Bolshoi catalog, we 
find 36,000 simulated halos that have exactly two satellites 
with Vmax > 50km i"' . These A^subs = 2 systems represent a 
first attempt at finding simulated halos that are analogs of the 
MW in terms of massive satellite content. What other proper- 
ties of the LMC and SMC might provide information on the 
mass and assembly history of the MW system? Repeated ob- 
servations over many years with HST and ground-based spec- 
troscopy have given us excellent limits on both the 3D posi- 
tion and velocity of both objects. Indeed, both of these prop- 
erties are significantly better constrained than are the circular 
velocities of these objects. We select simulated objects to be 
MC analogs based on Vmax, the maximum circular velocity of 
the object, ro, the distance of the object to the center of the 
MW, and s, the total speed of the object relative to the MW. 
These are summarized in Table [T] In order to be conservative 
with the uncertainties to account for systematics, we multiply 
the published formal errors in the speed by a factor of two 
when looking for MC analogs in Bolshoi (included in the er- 
rors shown in Table [1). This increase is n ecessary to bring 
the velo city measurements of t he SMC bv iKallivavalil et al] 
( I2006bl) and lPiatek et al] dlOOSh into agreement. 

Note that there is a slight discrepancy between the obser- 
vational measurements of satellite v^ax and the simulations: 
while observations measure the circular velocity curve for 
all components of the halo — dark matter and baryons — 
the Bolshoi simulation mod els only the dark matter content. 
iTruiillo-Gomez et"al] (120101) showed that ignoring the bary- 
onic component may cause a ~ 10% over-estimate in the mea- 
surement of Vmax for MC-sizcd objects. Since this correction 
is smaller than the error bars on the observations, we choose 
to ignore it. Additionally, there is some disagreement in the 
literature as to the allowable upper limit for Vmax- In partic- 
ular, some estimate s for the LMC are in excess of 100 km/s 
dPiatek et al.l 120081) . well above our adopted 80km/s Icr up- 
per limit. However, due to the rapi dly decUning abund ance 
of such massive subhalos (see, e.g.. iBusha et al.ll2010l) . our 
analysis is highly insensitive to changes in the upper error bar 
on Vmax- Finally, because of resolution effects, there is a ra- 
dial dependence to the incompleteness of satellite halos, with 
the large density contrast and limited force resolution mak- 
ing it difficult to identify subhalos closer to the center of their 
hosts. However, this incompleteness does not appear to corre- 
late with host mass, so we do not expect it to bias our results. 

4. INFERENCE 
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Fig. 1 . — Left: The MW mass inferred from the properties of its two most luminous satellites, the Magellanic Clouds. Lines show posterior PDFs (weighted 
histograms of Bolshoi halos) given information about (a) the existence of exactly two satellites with v^ax > 50 km s"' (blue), (b) the maximum circular velocities 
Vraax of the two satelHtes (red), (c) the distance of each satellite from the center of the MW (orange), (d) the speed of each satellite (green), and (e) all of these 
properties simultaneously (black). The combined properties define a sample of "satellite analogs" and give Mmw = 1-2^q 4 (stat.);';jj 3 (sys.) X IO'^Mq (68% 
confidence). The dotted fine shows a lognormal fit to this distribution, with parameters logmMMW = 12.2 ±0.1, and the dashed gray fines show the effect of 
bootstrap resampling: the uncertainty on the mean of the distribution from "sampling noise' in the catalogs was found in this way to be just 0.03 dex. 
Right: Comparison with various estimates for the mass of the MW from the literature. Dashed fines show results from: (a) the radial velocity dispersion profile 
(Battaglia et al. 2003[) (orang e); (b) the escape v elocity from halo stars |Smith et al. 20Q3) (red); (c) SDSS blue horizontal branch stars iXue et al. 2008) (green), 
and (d) the timing argument ILi & Whitell20Q8l) (blue). We assume lognormal error distributions, with asymmetric errors given by the quoted upper and lower 
confidence limits. The solid (black) line shows the posterior PDF for the MW mass from our satelfite analogs, and the dotted fine again shows thee lognormal fit 
to this distribution. 



TABLE I 

Observed properties of the LMC and SMC. 





LMC 


SMC 


Reference 


Vraax [km/s] 


65 ±15 


60± 15 


vdM02, S04, HZ06 


ro [kpc] 


50±2 


60±2 


vdM02 


i [km/s]* 


378 ± 36 


301 ± 104 


K06 



Note. — For a given satellite, Vraax is its estimated 
maximium circular velocity, its estimated distance from 
the Galactic center, and i its estimated speed relative to the 
Galact ic center. References are: vdM02 = van der Marel et aj] 
<2Q0g) ; S04 = [Stanimirovic et aU t2004D ; 

K06 = Kallivavi 3iretan i l2006all2006cl) 

HZ06 = IHarris & Zar itskv 1 20(^^ 

* En'ors on 1 have been increased relative to the published 
values for conservatism (see text). 

The Bolshoi halo catalog is large enough for some of its 
members to resemble the MiUcy Way system, with regards to 
its two most massive satellites. By weighting the Bolshoi ha- 
los according to how closely their satellites' properties match 
those of the MW's, we can infer the mass and assembly his- 
tory of the MW system. In this section we explain how this 
works, and derive the required weights from probability the- 
ory. 

4.L Observational constraints and the likelihood function 

As outlined above, the Bolshoi halo catalog can be thought 
of as a set of samples drawn from an underlying probability 
distribution. Each halo is characterized by a set of m param- 
eters, X, which includes the total mass of the halo and the 
properties of its subhalos, such as their masses, positions, and 

velocities: x = (Mvir,{v'max,'-o,i,...}LMC,SMc)- Given no ob- 
servational data, the set of Bolshoi halos provides a reason- 
able characterization of our prior PDF for the parameters of 
the Milky Way system, Pr(x). For example, if we were asked 
to guess the mass of the MW halo, we would do much worse 



by drawing a random number from some wide range, than we 
would do by drawing one Bolshoi halo at random from the 
catalog. 

However, we do have some observational data: we would 
therefore like to know the posterior PDF for the parameters 
of the Milky Way system, given this data d. Bayes' theorem 
shows this to be: 

Pr(x|d) cx Pr(d|x)Pr(x) (1) 

The first term on the right hand side is the joint likelihood 
function for the parameters, written as a function of the data. 
Given a particular parameter vector x, we can compute the 
value of this likelihood in the usual way, given assump- 
tions about the error distributions of the observational data. 
Our data d and their uncertainties a are given in Table [1] 
d= (v°nax,'"o''%'S°'") (for each MC). We use the superscript 
"obs" to make clear the distinction between the (constant) 
measured values, and the corresponding (variable) parame- 
ters of the model, which are the properties of the Bolshoi ha- 
los. We interpret the errors on the observational data as being 
Gaussian-distributed, such that for each datum di its likeli- 
hood function Pr((i, |x) is a Gaussian, with mean and standard 
deviation listed in Table [T] 

Since the = 6 observations were made independently by 
different groups using different techniques, the joint likeli- 
hood is just a simple product: 

N 

Pr(d|x) = nPrW|x), (2) 

= n ^(^max|Vmax,a2) (3) 
LMC, SMC 

y.N{rfyo.cTl)y.N{s°'''\s,cTl), (4) 
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where 



N{y\ti,cj) = 



1 



exp 



2a2 



(5) 



is the Gaussian probabiUty density in y with mean ^ and vari- 
ance (j^. 

Conditional on the values of the parameters, the data mea- 
surements are independent, and thus the likelihood function 
factors into the product of N=6 terms. In the prior, the dy- 
namical parameters and halo properties themselves are corre- 
lated, and not independent, because of the simulation physics. 
We have a measurement of satellite distance, r^^, and an inde- 
pendent measurement of satellite velocity s°'°^: the probability 
calculus makes it very clear that the PDFs for these two ob- 
served quantities can be multiplied together, even though the 
corresponding underlying parameters ro and s are not inde- 
pendently distributed. (We discuss correlations between the 
measurements of ro and s below.) The dynamics of satellites 
around halos is such that their positions and velocities are cor- 
related according to the orbits they are on. This is a good 
thing: it means that by measuring one we can learn about the 
other! In fact, it is this very interdependence that will allow us 
to infer the MW halo mass given measurements of its satellite 
properties. All the correlations between the model parameters 
are correctly taken into account by drawing samples from the 
joint prior PDF - that is, by using the Bolshoi catalog. 

4.2. Importance Sampling 

With the likelihood function in hand, we can now calcu- 
late the posterior PDF for the MW system parameters given 
our data. We have the prior PDF Pr(x) in the form of a set 
of n samples drawn from it; this is actually a very convenient 
characterization, and will allow us to derive an equally con- 
venient characterization of the posterior We compute poste- 
rior estimates using importance sar npling (see MacKav]| 2003l 
for an introducti on, and papers bv lLewis & Bridle.,20021 and 
iSuvu et ani20IOl for example applications in astronomy). In 
general, this technique involves generating samples from an 
importance sampling function, which are then weighted when 
computing integrals over the target PDF. In this paper, we 
choose the importance sampling function to be the prior PDF, 
so that the importance weights are proportional to the likeli- 
hood; this allows us to compute integrals over the posterior 
PDF, as shown below. 

These integrals include mean parameter values, confidence 
intervals and so on; a histogram of sample parameter values 
is a representation of the marginalized distribution for that pa- 
rameter, and is also a set of integrals (counts of samples in 
bins). 

In general, integrals over the posterior PDF can be written 
as follows: 



/(x)Pr(x|d)£/'"x: 



/ /(x)Pr(d|x)Pr(x)c/"'x 
/Pr(d|x)Pr(x)t/"'x 

E';/(X;)Pr(d|x,.) 
E';Pr(d|x,) 



(6) 



(7) 



where in the first line we have substituted for the posterior 
PDF using equation[T] and in the second line we have approx- 
imated the integrals with sums over samples drawn from the 
prior PDF. The denominator in each case is a normalization 
constant, /(x) is any function of interest to be integrated: for 
example, /(x) = x gives the posterior mean parameters. The 



highest importance samples correspond to halos that most re- 
semble, in terms of its bright satellite properties, the halo of 
the MW. 

4.3. Assumptions 

In this subsection, we examine the assumptions we have 
made in more detail. 

4.3.1. Data covariance 

As noted above, the constraints on ro and s are nearly, but 
not quite, independent. Working in galactic coordinates an d 
following the method outlined in Ivan der Marel et all JIOOl . 
the separation between the MCs and the center of the MW, ro, 
is given by the relation 



'•o = r0" + r; + 2(r0-r,), 



(8) 



where r0 is the location of the center of the galaxy and re 
location of the MCs, both relative to the sun. Similarly, the 
speed, s= |s|, for the MCs is given by 



s = v© + {^inVLn + ^iw^aw)re + v/„sU/„ 



(9) 



where v© is the relative motion between the sun and the center 
of the galaxy, ji^ and jiw are the measured north and west 
transverse components of the proper motion of the MCs on 
the sky, vios the measured line of sight velocity of the MCs, 
and un,Uw, U|os are the unit vectors defining the (north and 
west) transverse and radial directions of motion of the MCs 
relative to the sun in galactic coordinates. 

As can be seen in equations [8] and |9l r^,, the distance from 
the Sun to the MCs, is a necessary measurement for deter- 
mining both ro and s, hence these measurements are not fully 
independent. In order to determine the degree of dependance, 
we calculate COV(ro,s), the covariance between these vari- 
ables, and and show that this is significantly smaller than the 
measurement errors on the properties, 

COV(ro,5)<cr,„(T,. (10) 

The covariance can be explicitly written as 

COV(ro,.)=^;^a,,= 
are "fe 

1 r ■ r 1 

{re+^^)X —^ifll, + fll)X(Jr, (11) 

ro re 2y/s 



Using values reported in Ivan der Marel et all (120021) and 
lKallivavalil"etalT ( l2006al) . we measure COV(ro, s) = 6.52 and 
6.86 for the LMC and SMC, respectively. In both cases, this 
at least an order of magnitude smaller than the product ctj u ctj = 
72, and 208, allowing us to treat the measurements of ro and 
s as independent. 

4.3.2. Simulation time resolution 

The observational constraints on the positions of the MCs 
are actually tighter by a factor of two than the resolution 
of our simulations. The available Bolshoi outputs are such 
that a satellite with the L MC's radial velocity of ~ 90 km/s 
jKallivavalil et al.ll2006cl) will travel roughly 4 kpc between 
sequential snapshots, making this the effective uncertainty on 
the radial position of the simulated halos. We increase the size 
of the positional uncertainties accordingly, when calculating 
the likelihood function for the MC positions. 
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4.3.3. Importance sampling failure modes 

Typically there are two ways in which importance sampling 
can fail. The first is that the sampling function does not cover 
the domain of the target PDF, leaving parts of the PDF unsam- 
pled. In our case, we assert that the Bolshoi simulation does 
sample the parameter space, in that it contains (in large num- 
bers) halos and subhalos with masses, positions and velocities 
sufficiently close to the members of the MW system to make 
inferences meaningful. The second failure mode is that too 
few samples are drawn in the high importance volume, lead- 
ing to estimates dominated by a small number of high impor- 
tance samples. As we shall see, this "sampling noise" does 
lead to additional uncertainty on our inferences. Because of 
the very tight observational constraints on the properties of the 
MCs, relatively few Bolshoi halos receive significant impor- 
tance - we can only partially compensate for this by searching 
through multiple simulation outputs to improve the statistics 
of our sample. In Sections|5]and|6]below, we estimate the un- 
certainty, due to sampling noise, on each of our estimates by 
bootstrap resampling. 

4.3.4. The effect of the Galactic disk 

Finally, the lack of treatment of baryons in Bolshoi intro- 
duces a systematic error The more concentrated mass in the 
stellar disk at the halo center should increase the speed of the 
satellites orbiting their hosts at small radii. We can estimate 
the impact of this by artificially increasing the total veloc- 
ity of our simulated satellites by the circular velocity due to 
the stellar disk of the Milky Way, v^rc = ViGM^,/rsai), where 
M* = 6 X 10'°M^ (lKlvpin et al.ll2052l) . This correction is ad- 
mittedly simplistic, but is also small. To conservatively quan- 
tify the residual systematic error associated with this correc- 
tion, we take the difference in mass estimates with and with- 
out the above correction. This gives us an approximate up- 
per limit on the size of the two-sided systematic error due to 
baryon physics. 

5. THE MASS DISTRIBUTION OF THE MILKY WAY 

Figure [T] (left panel) shows the posterior PDF for the MW 
halo mass, computed by weighting every Bolshoi halo by the 
probability that its bright satellite population looks like that of 
the MW in various ways. Some observations provide signif- 
icantly more information about the host halo mass than oth- 
ers: the distance to and speed of the MCs are most constrain- 
ing, while the maximum circular velocity of the LMC and 
SMC provide almost no information (largely due to their mea- 
surement errors). Note that it is the combination of datasets 
that is most important, as internal degeneracies are broken, 
i.e., there is a high degree of covariance between positional 
and kinematic properties. The combination of all datasets 
gives Mmw = 1.2!|];^(stat.)!JS:^(sys.) x IO'^Mq (68% confi- 
dence) and a virial radius r,.,-,- = 250;!;3okpc. The systematic 
errors are estimated by repeating this process with and with- 
out the baryon correction discussed at the end of Section 4. 

How many halos contribute to this inference? We find 
114 1(7 matches (systems that are within an average of Icr 
from observations of the MCs in the 6 properties listed in Ta- 
ble [U and nearly 400 2a matches in the 60 snapshots. How- 
ever, many more contribute statistically. The effective num- 
ber of halos contributing to the posterior can be estimated as 
A/gff = A^/(l + Var(w)) (Neal 200Jj), where w is the normalized 
importance of each halo and A^, in our case, is 1.71 x 10^ (total 
number of halos overall 60 snapshots). We findA'eff = 10,051. 



While this is a healthy number of samples to compute statis- 
tics with, the low fraction of halos that contribute is driven 
largely by the combination of the total speed, s, and radial 
position, ro, constraints. This gives additional support to the 
idea that the MCs are currently at-or-near pericentric passage 
on their orbits around the Milky Way, since any object on an 
elliptical orbit spends such a small amount of time near peri- 
center Figure[T]shows the effects of the sampling noise on the 
inference of the MW mass: the faint blue-grey curves show 25 
bootstrap-resampled PDFs. The additional uncertainty on the 
central value is 0.03 in logjQ(MMw) - this is included this in 
the quoted error bars above, which were estimated from the 
sum of the resampled PDFs. 

Figure [T] (right panel) compares our result with previous 
MW halo mass estimates from the literature. The LMC and 
SMC properties lead to a halo mass that is in excellent agree- 
ment with the dynamical estimates in the literature. In partic- 
ular, our results are in near perfect ag reemen t with the most 
recent stellar velocity measurements (IXue et al. 2008), with 
similar error bars. 

Throughout the paper, we refer to the collection of hosts 
weighted by the likelihood of the v,nax, A), and s of their satel- 
lites as "satellite analogs" of the MW. Note that this does not 
imply that we have selected a specific subset of hosts. Rather, 
we have taken all hosts with Nsuhs = 2 and weighted each ob- 
ject by its satellite properties. It is this sample of weighted 
objects that defines our satellite analogs. It is worth noting, 
however, while we formally use all 35,000 Nsubs = 2 halos in 
Bolshoi, the ^ 500 1- and 2-a halos provide more than 95% 
of the total weight. 

Because we want to further understand what impact the 
MCs have on other other halo properties, we also define "mass 
analogs" to be the set of halos randomly drawn from the mass 
PDF of the satellite analogs. Thus, the mass analogs have the 
same PDF of virial masses as the black line of Figure[Tl but no 
constraints on their satellite properties. Comparison between 
our satellite analogs and mass analogs allows us to disentan- 
gle impacts on the system due to the satellite properties from 
those due to the particular mass range probed by our satellite 
analogs. 

When interpreting Figure [T] it is also helpful to consider 
the full dependence of the predicted satellite properties on the 
mass of the host halos. This is shown in Figure|2l which com- 
pares the satellite parameters v,„oi , ro, and s with the Myir of 
their hosts. The plots show contours containing 68, 90, 95, 98, 
and 99% of our prior probability (that is, all hosts with exactly 
two satellites, black lines), as well as the exact locations of our 
2-(T (open red circles) and l-cr hosts (filled circles). Also plot- 
ted are the observed properties of the LMC and SMC (blue 
and green lines). These plots highlight a number of important 
trends. First, we can see that the MCs are atypical subhalos in 
most regards. The LMC in particular is roughly a 2-a outlier 
in each of these properties. Second, correlations between pa- 
rameters can be seen, such as the that between speed and host 
mass, which shows explicitly how the high observed speed 
of the satellites skews the posterior PDF towards higher mass 
hosts. Additionally, it is interesting to note that, while rare, 
there are a few objects with Myir < 1O"M0 with satellites that 
are well matched to the MCs. These satellites are fast mov- 
ing, massive subhalos roughly 50 and 60 kpc from their host 
centers - well within their host virial radii, but not energeti- 
cally bound to their hosts, and hence merely transient events. 
Finally, while it can be seen that observations of s provide the 
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Fig. 2. — Relationships between Myirhost and various satellite parameters: 
Vmax (top left), Co (top right), and s (bottom left). In all panels, the black 
contours denote the regions containing 68, 90, 95, 98, and 99% of our prior 
distribution, all satellites in a host with W,,,;, = 2. The open red and filled 
circles denote the satellites for 2- and l-tr halos, those that contribute most to 
the posterior PDF. Blue and green lines show the observed values and ± 1 - cr 
uncertainty for each property, for the LMC and SMC respectively. Similarly, 
the blue and green circles denote the values for our identified LMC and SMC 
analogs. 



most Stringent constraints individually, these plots further em- 
phasize that it is the combination of observed properties that 



is necessary to place tight constraints on Mmw- 

We can also determine the impact of the MCs on the in- 
ternal mass distribution of a halo by comparing the density 
profiles for our satellite and mass analogs. This gives us a 
handle on how typical the MW is for a halo of its mass. We 
find that the presence of the MCs has only a modest impact on 
the halo concentration, but can impact density profiles signif- 
icantly in other ways. In particular, when looking at the mass 
enclosed within a fixed 8 kpc (the distance from the sun to 
the center of the galaxy), satellite analogs have a 60% higher 
central density than mass analogs. This tells us that the MCs 
are correlated with a more strongly peaked inner dark mat- 
ter distribution. For the full density profile (which includes 
contributions from substructures), however, satellite analogs 
have concentration c,,„- = 11 ±2, a little less than l-cr higher 
than c,,,> = 8.7 ± 3.5 for the mass analogs (here, Cviv- = '"n>/'".s)- 
While this still implies a correlation between the presence of 
the MCs and a more peaked mass distribution, it shows that 
the impact is much weaker at larger radii than it is for the 
central core. 

6. THE ASSEMBLY OF THE MILKY WAY HALO 

We now turn our attention to the assembly history of the 
Milky Way. We use the same importance-sampled halos from 
the previous section to infer the MW assembly history, in the 
same way as we inferred its mass and density profile. Fig- 
ure[3]shows the distribution of accretion times for these satel- 
lite analogs (where tacc is the time since the subhalo first came 
within 300 kpc of the host); we also show the accretion time 
PDF for all hosts with A^subs = 2, and for the mass analog sys- 
tems defined above. The mass analogs (red line) clearly show 
two populations. The first consists of halos whose subhalos 
were accreted at high redshift, when the host halo was in its 
exponential g rowth phase, which suppresses tidal stripp ing of 
the subhalos dWechsler et al.ll2002l: iBusha et al.ll2007l) . The 
second population consists of halos with recently accreted ob- 
jects that have not had enough time to undergo significant tidal 
disruption. The relative size of these populations changes 
when we apply the observational likelihoods. Requiring that 
a host has A^sub,s = 2 massive subhalos (with unconstrained 
speeds and distances, dotted line) has little impact. However, 
for satellite analogs (black line), the size of the recently ac- 
creted population increases dramatically. This is primarily 
driven by the combined requirement that the satellites have 
both a high radial velocity and are close to the center of the 
halo, and argues that there is roughly a 72% chance that MCs 
are recent arrivals, accreted within the past 1 Gyr While it 
is clear that low number statistics are strongly impacting the 
shape of the probability distribution for the accretion time of 
the mass analogs in Figure [3] our bootstrap analysis shows 
that the measurement of 72% of all subhalos accreting within 
the last 1 Gyris a robust result that is not sensitive to the small 
number statistics. This is again demonstrated by the gray -blue 
lines in Figure [3] which show the distriubtion of accretion 
times in 25 or our resamplings. In all of these cases, there is a 
strong peak at recent accretion times. Our bootstrap analysis 
shows that the measurement of 72% of all subhalos accreting 
within the last 1 Gyris a robust result, that is not sensitive to 
sampling noise. Indeed, in 95% of our resamples, there is a 
> 65% probability that the Magellanic clouds were accreted 
within the past 1 Gyr. 

Did the MCs arrive together? Figure|4]shows the difference 
in accretion times, Afacc, for the two most massive satellites of 
all hosts with A^subs = 2 and for the two satellites of the satellite 
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analog hosts. For satellite analogs, the distribution is strongly 
peaked towards small Afacc, with roughly 50% of satellites 
having been accreted simultaneously (within a Gyr of each 
other). These simultaneous accretions correspond very tightly 
with the recently accreted population. The noise in this plot, 
and in particular the peak around At^^c ~ 8 - 9Gyr is driven 
by a few very well matched (high weight) halos that had one 
satellite accrete within the last Gyr and the other early on dur- 
ing the exponential buildup phase (see the weaker secondary 
peak in Figure O. This highlights the current level of noise in 
our analysis due to our modest statistical size. We anticipate 
that, with better statistics, this high-Afacc bump will smooth 
out, making a smoother distribution that is even more sharply 
peaked at Afacc = 0. Still, using out bootstrap analysis, the 
large fraction of objects accreting within a Gyr of each other 
appears to be quite robust, as can be seen by the locus of blue- 
gray lines in Figure |4] For comparison, halos with A^subs = 2, 
shown as the dashed line, have a much weaker preference for 
simultaneous accretion. 

Both this result and that of Figure [3] favor a method for the 
creation of the Magellanic Stream other t han tidal di s ruptio n 
by the MW. For example, the model of iBesla et al.l (120101) . 
who found good agreement with the dynamics of the Stream 
for a model in which it was created by tidal disruption of 
the SMC by the LMC before they were accreted as a bound 
pair Additionally, the satellites in the satellite-analogs have a 
high degree of spatial correlation, with a mean separation of 
48 ± 8kpc at the epoch where the simulations best match ob- 
servations of the MCs, about 3a larger than the observed MC 
separation of 25 kpc ( iKallivavaUl et al.ll2006ah . yet roughly 
3(7 closer than expected for two randomly drawn points at the 
appropriate radii for the MCs. This provides some support for 
the notion that the MCs may be a bound pair. However we 
do not detect a lower relative speed of the two subhalos in the 
satellite analog sample than in the mass analog sample. The 
Bolshoi simulation lacks the volume to address this question 
more directly: ideally, we would like to further weight our 
sample by the observed separation between the two satellites, 
but this would reduce the effective number of halos below the 
sampling noise limit. While we cannot say very much about 
the boundedness of the LMC and SMC pair, we can make the 
following point: that even without using the observed angular 
separation of the LMC and SMC as a constraint, we find a 
high probability of them arriving at the same time. 

7. CONCLUSIONS 

The advent of high-resolution cosmological simulations 
which sample the dynamical histories of large numbers of 
dark matter halos in a wide range of environments provides us 
with a new approach for determining the properties of individ- 
ual halos given their observational characteristics. Here, we 
use the observed properties of the Magellanic Clouds to con- 
strain the mass distribution and assembly history of the Milky 
Way. In comparison to previous efforts which use detailed 
observations but a necessarily simplified dynamical model, 
our approach uses simple observations and statistical infer- 
ence from sampling of a detailed and cosmologically consis- 
tent dynamical model. 

Our principal conclusions are: 

I. We infer the MW halo mass to be Mmw = 
L2!|];^(stat.)![!:^(sys.) x IO'^Mq (68% confidence), in 
very good agre ement with the recent stellar velocity 
measurements (IXue et al.l 120081) . with similarly sized 
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Fig. 3. — Posterior distribution of satellite accretion times, from the MW 
satellite analogs (black), from hosts with exactly two subhalos (dashed), 
and from MW mass analogs (red). Selecting hosts with MC-like satellites 
strongly weights the distribution towards recent accretion. The blue-gray dot- 
ted lines show the distributions for 25 random MW satellite analogs drawn 
from our bootstrap analysis. 

error bars. 

2. The MW halo has a slightly higher concentration than 
is typical for its mass: 1 1 ± 2, compared to 8.7 ± 3.5. 
Additionally, the density within 8 kpc is 60% higher for 
satellite analogs than for mass analogs. 

3. Typical bright satellites of halos with MW mass were 
accreted at a range of epochs, generally at high redshift 
(c. 10 Gyr ago) or much more recently (within the last 
2 Gyr). Because of their high speed near the center of 
the halo, we find a 72% probability that the Magellanic 
Clouds were accreted within the last Gyr. We also find 
a 50% probability that the MCs were accreted within 
1 Gyr of each other 

This approach clearly allows one to explore a wide range 
of additional properties of MW-like halos; however, it places 
challenging requirements on the simulations used. In particu- 
lar, for MW studies we are still limited by the relatively small 
simulation volume used. With the constraints used here, we 
found just one good fit halo per ~ 500Mpc^, emphasizing the 
large volume required to perform this analysis. As more crite- 
ria are applied we will need to sample a larger range of forma- 
tion histories and environments drawn from a larger cosmo- 
logical volume. In addition, the properties of its smaller satel- 
lites are not accessible with our present resolution. Pushing 
forward on both simulated resolution and volume will be es- 
sential to realize the full potential of this approach. Addition- 
ally, the simulations used here ignore the impact of baryons on 
the dark matter distribution. We have included a simple model 
of the stellar disk which is applied to the satellite velocities to 
get a handle on the impact of this systematic, but the model is 
simplistic and the effects of the baryonic component need to 
be further studied. 

As we were completing this work, iBovlan-Kolchin et al.l 
( 120101) presented results from a similar study. Their princi- 
pal results regarding the mass of the MW and the accretion 
history of the LMC and SMC are in reasonable agreement 
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Fig. 4. — Difference between accretion times for the two most massive 
satellites in MW-like systems. The solid line represents satellite analogs; 
the dashed line shows all systems with exactly two subhalos. The blue-gray 
dotted lines show the distributions for 25 random MW satellite analogs drawn 
from our bootstrap analysis. 

with our own, although they favor a somewhat larger MW 
mass, Mvii- ~ (2-3) x IO'^Mq. These studies differ in the 
16 times larger volume simulation used here, the cosmology 
of the simulations, as well a s our use of statistical inference. 
iBovlan-Kolchin et al.l (120101) used a simulation with cosmo- 
logical parameters taken from the 1-year WMAP results, with 
fijw = 0.25, a% = 0.9, and h = 0.73, while the Bolshoi sim- 
ulation uses the more recent WMAP7 results and features 
a l ower = 0-82, wi th 51,,, = 0.27 and h = 0.7. Adopting 
the [Tinker et"aLl ( l2008l) mass function, we can understand the 
impact of these cosmological differences by considering the 
abundances of MW and MC mass objects. Differences come 
from a number of competing effects. First, due to cosmologi- 
cal differences, their lower cts will tend to suppress our mass 
function, while the lower h will increase the halo masses. In 
this case, the differences in h has the more significant impact, 
and causes the abundances of halos withMyi, = 1.2 x lO'^M© 
to be suppressed by about 10% in the Millennium cosmology. 



However, the more relevant number is the ratio of the num- 
ber density for our selected satellites to their selected hosts. 
Adopting the fiducial mass for the LMC, 100 times lower 
than that of the MW, we see that the ratio in number densi- 
ties is 65.0 in Bolshoi vs. 63.2 in Millennium, a suppression 
of roughly ^ 3%. Because the mass function is relatively flat 
here, reproducing the Bolshoi abundance ratio requires that 
the host mass increases to 3.6 x 10'~Mq, a correction that 
brings our results into better agreement with theirs. However, 
it is also important to note that significantly different selection 
criteria for identifying "MW -like" objects were employed. 
iBoylan-Kolchin et al.l (1201 Ol) selected hosts whose two largest 
subhalos were within 0.75 rioo and had similar total velocities 
and stellar masses to the MCs using abundance matching to 
estimate the stellar content of their simulated halos. In this 
work, we select objects with exactly two subhalos more mas- 
sive than Vniax = 50km i"' within a fixed 300 kpc aperture and 
then weight our sample according to how well the subhalos 
look like the MCs in terms of v^ax, position, radial velocity, 
and total speed. These differences likely account for the ten- 
sion in the resulting Mmw PDFs. In particular, their focus on 
the total speed, s, of the magellanic clouds has likely biased 
their mass estimates for the MW high. This is shown by com- 
paring the green and black lines in Figure [T] which show the 
likelihood distribution of masses for hosts having two satel- 
lites with similar speed to the MCs (green), as well as for 
hosts with two satellites with similar speeds, masses, and po- 
sitions (black). The green di stribution presents a much be tter 
agreement with Figure 13 of lBovlan-Kolchin et al.l (120101) . 
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